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We study the quantum phases of mixtures of ultra-cold bosonic atoms held in an optical lat- 
tice that confines motion or hopping to one spatial dimension. The phases are found by using 
Tomonaga-Luttinger liquid theory as well as the numerical method of time evolving block decima- 
tion (TEBD). We consider a binary mixture with repulsive intra-species interactions, and either 
repulsive or attractive inter-species interaction. For a homogeneous system, we find paired- and 
counterflow-superfluid phases at different filling and hopping energies. We also predict parameter 
regions in which these types of superfiuid order coexist with charge density wave order. We show 
that the Tomonaga-Luttinger liquid theory and TEBD qualitatively agree on the location of the 
phase boundary to superfiuidity. We then describe how these phases are modified and can be de- 
tected when an additional harmonic trap is present. In particular, we show how experimentally 
measurable quantities, such as time-of-flight images and the structure factor, can be used to dis- 
tinguish the quantum phases. Finally, we suggest applying a Feshbach ramp to detect the paired 
superfiuid state, and a 7r/2 pulse followed by Bragg spectroscopy to detect the counterfiow superfiuid 
phase. 
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I. INTRODUCTION 

Bose-Einstein condensation [l| is a fascinating many- 
body phenomenon. It demonstrates the significance of 
quantum statistics at low temperature. Identical bosons 
can occupy the same single particle state and are in fact 
more likely to do so than classical particles. At a critical 
temperature, a gas of bosons undergoes a phase transi- 
tion towards a state in which a macroscopic fraction of 
the particles occupy the lowest energy state, creating a 
condensate. Such a state was realized in ultra-cold atom 
systems in f3| , demonstrating that the technology of cool- 
ing and manipulating atoms had reached a level of control 
with which novel states of matter could be generated and 
studied. 

In the case of a Fermi gas, the Pauli exclusion prin- 
ciple prevents such a phenomenon to occur, because 
no single particle state can be more than singly occu- 
pied. However, the phenomenon of condensation can 
still occur in Fermi systems via a different mechanism: 
fermions can form pairs to create composite bosons. The 
bosonic particles then form a condensate of pairs. Con- 
ventional superconductors, for example, were understood 
as a condensate of electron pairs . In ultra-cold atoms, 
fermionic condensates of this type were created in Q • 

Interestingly, this mechanism of condensation of pairs 
is not limited to fermionic systems but can occur in 
bosonic systems as well. In fermionic systems, formation 
of Bosonic pairs necessarily occurs before condensation. 
In bosonic systems this mechanism can be favored ener- 



getically, and will typically be in competition with single 
particle condensation. 

In two types of composite bosons were predicted 

for a binary Bose mixture in a optical lattice: pairs and 
anti-pairs. For attractive mutual interactions, a bosonic 
mixture can form pairs of atoms which then form a paired 
superfiuid (PSF) state, as is visuaHzed in Fig. [TJ For 
repulsive interactions, at special fillings, the atoms can 
form anti-pairs, which can be interpreted as pairs of one 
atom of one species and one hole of the other species. 
These anti-pairs can then generate a counterfiow super- 
fiuid (CFSF) state, visualized in Fig. H Most of their 
simulations were performed for two dimensional systems. 

Quantum phases of atoms in optical lattices have been 
experimentally studied. Following the prediction by 
Jaksch et al. in Q, the Mott insulator (MI) to super- 
fiuid (SF) transition was realized in Ref. Q in a three 
dimensional lattice. In p| this transition was achieved in 
ID. More recently, Ref. |lO| observed the two dimensional 
(2D) transition. 

In one-dimensional gases quantum phases have quasi- 
long range order (QLRO), rather than true long range 
order. QLRO of an operator 0{x) is defined as follows: 
The correlation function R{x) = {O^x)O{0)) falls off 
algebraically as R{x) ^ as \x\ oo with a > 0. 

Various order parameters 0{x) will be defined in the text. 
In contrast in higher dimensional bosonic systems corre- 
lation functions can have true long range order, where 
correlation functions approach a finite value. Power-law 
scaling in a ID optical lattice has been observed in 
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Figure 1: Sketch of a condensate of pairs. Atoms of each 
species (red/green) pair together and form a paired superfluid 
(PSF) state. 



Figure 2: Sketch of a condensate of anti-pairs. Here, atoms 
of one species are strongly anti-correlated with atoms of the 
other species, creating a counterflow superfluid (CFSF) state. 
These composite bosons can also be thought of as a pair of 
one atom of one species and one hole of the other species. 



They observed the Tonks-Girardeau regime of strongly 
interacting bosons. 

In this paper we consider a two-component Bose mix- 
ture held in an optical lattice that only allows atoms 
to hop in one spatial dimension. We ask the question 
of how the superfluid as well as other phases or orders 
can be reaHzed. We assume that the two species of 
the mixture have the same filling v, restricted to the 
range < < 1. The phase diagram of these mix- 
tures is determined using Tomonaga-Luttinger Hquid the- 
ory , which gives the universal phase diagram in terms 
of a few effective parameters. Based on the univerisal 
phase diagram, we generate the numerical phase dia- 
gram using the time-evolving block decimation (TEBD) 
method [IB, [13, [iBl • With these two approaches we 
find that CFSF can exist for 1/ = 1/2 (half-filling) and 
repulsive interaction, whereas PSF can exist for v < 1 
and attractive interaction (see also [l^). 

We also find that charge density wave (CDW) quasi- 
order can coexist with both PSF and CFSF, as well as 
single particle superfluidity (SF). The regimes in which 
CDW and SF quasi-order coexist constitute a quasi- 
supersolid phase (20l . [2l| . Similarly, the regimes where 
CDW and PSF quasi-order coexist is a quasi-supersolid 
of pairs and in the case of CFSF, a quasi-supersoHd of 
anti-pairs. Previous work has predicted coexistence of 
CDW and PSF for ID Bose mixtures [13, [3 and bi- 
layer 2D lattice bosons with long-range interactions [23 |. 
and that of CDW and CFSF for ID Bose-Fermi mix- 
tures [2i|,[2i|. 



We then address the question whether PSF and CFSF 
can be realized and detected in experiment. To simulate 
the effect of a global trap, we numerically study a mixture 
confined by a harmonic trap and find that PSF and CFSF 
can indeed exist in such trapped systems. Their existence 
can be detected through various measurements. The PSF 
phase can be detected by using a Feshbach ramp, simi- 
lar to what has been used in BEC-BCS experiments 01, 
which generates a quasi-condensate signal in the resulting 
molecules. The CFSF phase can be detected by applying 
a 7r/2 pulse followed by Bragg spectroscopy. This gen- 
erates a quasi-condensate signal in the structure factor. 
Time-of-fiight expansion can also be used to show the ab- 
sence of single particle superfluidity in PSF and CFSF. 
Measuring the structure factor via Bragg spectroscopy 
can be one way of detecting CDW order. 

This paper is organized as follows: in Section[lll we in- 
troduce the model that is used to describe the system; in 
Section llllj we use Tomonaga-Luttinger liquid theory to 
derive the phase diagram. The numerical approach and 
results are discussed in Section IIVI Speciflcally, phase 
diagrams of the homogeneous system are presented in 
Section IIV A| and the realization and detection of PSF 
and CFSF are discussed in Section HV Bl We conclude in 
Section |Vl 



II. HAMILTONIAN 

Ultra-cold bosonic atoms in optical lattices can be well 
described by Bose Hubbard models Here, we con- 
sider a mixture of two types of atoms confined to a one- 
dimensional lattice system. The Hamiltonian of such a 
system is given by 

JV-l N 

H = ~t {bi,iba,i+i + h.c.) + U12 ni^in2,i 

a=l,2 i=l 



u 



N 



-T ^ ^?^a,^(^^a,i - !)■ 



(1) 



a=l,2 i=l 



We denote the different types of atoms with index a = 
1,2, and the lattice site with index i. We assume that the 
two species have equal particle density v < \^ the same 
intra-species interaction [/ > and hopping parameter 
t > 0. The inter-species interaction is given by Ui2. The 
operators b\ ^ and ba.i are the creation and annihilation 

operators for atoms of type a and site i and Ua,: ~ 
are the number operators. 
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III. TOMONAGA-LUTTINGER LIQUID 
APPROACH 



The universal behavior of this system can be found 
within a Tomonaga-Luttinger liquid description In 
this paper, we are interested in the phase diagram of the 
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system at various densities and interactions. First, we 
switch to a continuum description, 6a,i ba{x), and ex- 
press the operators bgjx) through a bosonization identity, 
according to Haldane [25l . [2^ : 



ba{x) 



(2) 



where the real-space density of each species is n = ly / ql 
and ai is the lattice constant. The lattice sites are at 
positions x = iaL- This expression is a phase-density 
representation of the Bose operators, in which the square 
root of the density operator has been written in an in- 
tricate way. The fields IIi.2{x) describe the small am- 
plitude and the long wave length density fluctuations. 
The fields 61,2 (a^) are given by 61,2(2;) = nnx + 6*1,2(2;), 
where ^1,2(2;) — tt dylli^2{y)- The fields (l)i,2{x) de- 
scribe the phase, and are conjugate to the density fiuc- 
tuations ni,2(a;). 

The contact interactions between the densities in [T] 
written in Haldane's representation generate an infi- 
nite series of terms that contain exp(2mii(7rna; -I- 0i) + 
2m2i(T^nx + 62)) J where mi and m2 are some integers. A 
term of this form can only drive a phase transition, if the 
oscillatory part 27rminx-t-27rm2na; vanishes for all lattice 
sites. This leads to the requirement miv -\- m2V = ms, 
with 7713 another integer [23|. As a further requirement, 
small integers mi and m2 are necessary, because the scal- 
ing dimension of the term scales quadratically in mi and 
m2. 

For the range < < 1, we find that there are three 
different cases: unit-filling {v = 1), half-filling {v = 1/2), 
and non-commensurate filling (v ^ 1 and ^ 1/2). It 
can be checked, using renormalization group arguments 
as below, that higher forms of commensurability do not 
generate new phases, but that either phase separation 
or collapse is reached first. Our numerical findings are 
consistent with this. 

Non- commensurate filling. The action of the sys- 
tem, assuming a short-range spatial cut-off ro, at non- 
commensurate filling is given by [l^, [l^, [2^ : 
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Ui2aL 



dxdld: 



•x02 



(2^ro)2 



cos(26li - 26*2)] (3) 



The first line of the action is characterized by a Luttinger 
parameter K and a velocity v, contained in r = {vt,x). 
This part of the action, without the coupling between the 
two fields Oa{x), generates a linear dispersion w = v\k\, 
where . v should therefore be interpreted as the phonon 
velocity. The Luttinger parameter if is a measure of the 
intra-species interaction U. We will be interested in the 
regime U >t, in which we have approximately [23| 



K 



1 + 



8t sin irv 



U 



(4) 



The velocity v can also be related to the parameters of 
the underlying Hubbard model by 



V « Vp{l — 8tv COSTTV /U) 



(5) 



where is the 'Fermi velocity' of an identical system 
of fermions, vp = 2{aLt/h)svin:v^ and kp is the 'Fermi 
wave vector', kp — irn. Here, h is the Planck constant. 

The two fields Oa{x) are coupled by the inter-species 
interaction. The interaction term Ui2nin2 in the under- 
lying Hubbard model generates both the term containing 
dxOidx02, as well as the backscattering term (l3l.[2^ con- 
taining cos(26'i — 2^2)- The action S is only well-defined 
with a short-range cut-off tq. It is proportional to 1/n. 
At this scale, ga is approximately given by 



ga = Ui2aL/{vh). 



(6) 



We diagonalize the quadratic part of the action by 
switching to the symmetric and antisymmetric combi- 
= -75 (6*1 ±6*2)- For the two sectors we find 



nations 



'S/A 



K 



S/A 



{l/K^ ±Ui2aL/{vhTTK)) 



-1/2 



(7) 



as effective Luttinger parameters. To lowest order in U12 
this gives Kg/A ~ K ^ Ui2aLK^ / (2TTvh). The effective 
velocities are vs/a = vy/l ± Ui2aLK/{nvh). Collapse 
(phase separation) of the superfiuid phase is when vg/A 
is imaginary. We note that Ks diverges when collapse 
(CL) is approached, and that Ka diverges as the system 
approaches phase separation (PS). 

The anti-symmetric sector contains the nonlinear 
backscattering term cos{2\/29a)- To study its effect, we 
use an RG approach. We renormalize the short-range 
cut-off ro to a slightly larger value, and correct for it at 
one-loop order. The resulting fiow equations are given 
by [H: 



^ = {2-2KA)g. 



dKA 

dl 



9c, 7^3 



The fiow parameter I is given by 

I = loge f ^ 



(8) 
(9) 

(10) 



where Tq is the new cut-off that has been created in the 
RG process. 

The fiow equations [8] and [9] have two qualitatively dif- 
ferent fixed points: Either g^ diverges, which in turn 
renormalizes Ka to zero, or g^ is renormalized to zero 
for finite Ka = K\. In the latter case, the action S is 
quadratic in 63 and 9a. For the parameter Ks, we use 
the bare value given in Eq. [7l 

As mentioned in the introduction, we can determine 
the phase diagram by studying the long-range scaling 
behavior of correlation functions, {0\x)0{y)) , of var- 
ious order parameters 0{x). In particular, the single- 
particle superfiuid order parameter is Osp = ba{x) with 
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Figure 3: Phase diagram of a bosonic mixture at non-unit and 
non-half-filling. For attractive interactions U12 and K < 2 the 
system can form a paired superfiuid state, in the regime la- 
beled PSF and PSF(CDW). This phase can coexist with CDW 
order for weaker interactions. For large repulsive (attractive) 
interactions f/12 the system phase separates (PS) (collapses 
(CL)). For the remaining regime the system shows single par- 
ticle superfiuidity (SF). This can coexist with CDW order, 
resulting in a quasi-supersolid (SS) regime. 



a = 1,2. The CDW order is related to the 2kF wavevec- 
tor component of the density operator, Ocdw — n-a. 
PSF is described by Opsf = bi{x)b2{x), and CFSF by 
OcFSF = h\{x)h2{x). In the homogeneous system, it 
suffices to study 

G{x) = {hl{x)bam,a^l,2 (11) 

Rn,a{x) = (na(x)na(0)),a= 1,2 (12) 

Rs{x) = (&l(x)4(a;)6i (0)62(0)) (13) 

Ra{x) = {b\{x)b2{x)hmlm. (14) 

We find that away from collapse (CL) and phase sepa- 
ration (PS), the correlation functions scale either alge- 
braically or exponentially. For algebraic scaling, we have 



G{x) 

Rs{x) 
Ra{x) 



asF = 2 - l/{4Ks) - 1/{AKa) (15) 

cos(2fcFa;)|a;|"'=^"^^"^ 

acDW = 2~Ks-Ka (16) 

\x\"^'^-^,apsF = 2-l/Ks (17) 

|^|acrs^-2^ = 2 - 1/Ka. (18) 



where the scaling exponents ao are determined by Ks 
and Ka after the RG fiow. For the case that diverges 
in Eqs. [8] and [9] and Ka is undefined, these expressions 
can still be used. We set Ka to zero, and find that acDW 
and apsF are well defined. Hence Rn,a and Rs still show 
algebraic scaling. On the other hand, acFSF and asF 
become —00 and G and Ra scale exponentially. 

We can identify regimes where different scaling expo- 
nents are positive based on the relationship between the 



scaling exponents and Ks/a after the fiow. This de- 
termines the different quasi-long range orders that are 
present. The resulting phase diagram is shown in Fig. [31 
as a function K and Ui2aL/{vh), as appearing in the ac- 
tion in Eq.[3l These two parameters determine the initial 
values of the fiow equations through equations [7] and [H 
We can estimate the phase boundary between PSF and 
SF. For small Ui2aL/{vh) this boundary is near the point 
Ka ~ 1 and ~ 0. For that limit, Eq. ([9]) can be 
linearized to 



dKA 
dl 



9a 

27r2 



(19) 



and the expression A — tt'^{1 — Ka)^ — 5^/4 becomes 
an invariant of the fiow. From the properties of the RG 
fiow of a Berezinskii-Kosterlitz-Thouless transition (see 
e.g. [3, H^l), the phase boundary is given by ^ = and 
Qa < 0. Using the expressions of Ka and v in terms of the 
Hubbard parameters, we estimate the critical interaction 
U12 for PSF to occur at 



U12 
U 



= -32— ^sin2(7rzy) 



(20) 



The phase boundary between supersolid (SS) and SF has 
been derived in Ref. jlBl- 

Half -filling . In the case of half-filling, another non- 
linear term has to be introduced in the action 

Suk = J d\cos{29i+ 292). (21) 

This term describes Umklapp scattering. At the initial 
cut-off ro ~ 1/rt, guk is approximately given by UuaL/v. 
In addition to the RG fiow in the antisymmetric sector 
we now also have 



dguk 
IF 
dKs 
dl 



= {2-2Ks)guk 



9lk 



2-k' 



2^S 



(22) 
(23) 



in the symmetric sector. Proceeding along the same lines 
as for the non-commensurate case, we find the phase di- 
agram shown in Fig. [H 

We estimate the SF-CFSF phase boundary in the same 
way as the PSF-SF boundary. We find 



U12 

U 



32— sin^(7riy) 



(24) 



Unit-filling. At unit-filling we have to introduce a term 
of the form 



Si 



2gi 



dV(cos(26'i) +cos(202)) • (25) 



(2^ro)2 

The resulting RG fiow for this system is given by 



dguk 
dl 



{2-2Ks)guk + a3 



gUKA-Ks) 
2n 



(26) 
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Figure 4: Phase diagram of a bosonic mixture at half-filling. 
In addition to the phases that appear in Fig. O the system 
now develops a counterfiow superfiuid (CFSF) phase, which 
can coexist with CDW order. 





Rs{x) 


Ra{x) 


Gix) 


MI 


Exp. 


Exp. 


Exp. 


SF 


Alg. 


Alg. 


Alg. 


CFSF 


Exp. 


Alg. 


Exp. 


PSF 


Alg. 


Exp. 


Exp. 


CL/PS 


Rs{x), Ra{x) undefined 



Table I: Definitions of Mott insulator (MI), superfiuid (SF), 
counterfiow superfiuid (CFSF) and paired superfiuid (PSF) 
orders in terms of the long-range behavior of the correlation 
functions Rs{x), RAix), and G{x) . Each of these can ei- 
ther show algebraic (Alg.) or exponential (Exp.) decay when 
the system is away from collapse (CL) or phase separation 
(PS). From the Tomonaga-Luttinger liquid theory, Rs{x) and 
Ra{x) approach a constant (or Ks/a diverges) when the sys- 
tem approaches CL/PS regime. For the numerical calculation 
in the CL/PS regimes, the behavior of the correlation func- 
tions is inconclusive and we assign the phase from additional 
observables as discussed in the text. 



dga_ 
dl 
dgi 
dl 

dKA 
dl 

dKs 
dl 



(2 - 2KA)ga + as 



gl(Ks-KA) 



Ks+Ka , 
(2 + ag 



27r 

gukKs + gcrKA 



(27) 
)3i(28) 



9uk_j^3 



{Ks + Ka)KI (29) 



27r2 



s 



167r' 



{Ks + Ka)KI 



(30) 



where as is some non-universal parameter [28l|. The be- 
havior of this set of equations depends strongly on the 
initial value of gi. For small values of gi, four phases 
can be stable: Single-particle superfluidity, CFSF, PSF 
and a Mott phase. For large values only single-particle 
SF and MI are stable. We determine with our numerical 
approach, that the Hubbard model falls into the second 
category, i.e. there is only a single-particle SF and a Mott 
state at unit-filling. 

Having established the universal behavior of the sys- 
tem from Tomonaga-Luttinger liquid theory, we now 
want to connect the phase diagram with the parame- 
ters in the Hubbard model. The expressions |4] and [H 
which relate the Luttinger parameter K and the velocity 
V to microscopic parameters of the Hubbard model, are 
only approximate, no full analytic expression is known. 
In addition, only some phase boundaries are predicted 
reliably, because we use perturbative RG in the gcr. We 
expect that the analytic calculation only predicts the gen- 
eral structure of the phase diagram, as well as the de- 
cay behavior of the correlation functions. To obtain the 
phase diagram in terms of the parameters in the Hubbard 
model, we need to use numerical methods. The next sec- 
tion describes the numerical determination of the phase 
diagram. 



IV. NUMERICAL APPROACH 

We use the time-evolving-block-decimation (TEBD) 
method to study our discrete one-dimensional two- 
species Hubbard Hamiltonian. With this method, ex- 
plained in Appendix[Al we obtain an approximate ground 
state solution. We consider N lattice sites with hard-wall 
boundary conditions and express the Hubbard parame- 
ters in units of the intra-species interaction U. The num- 
ber of sites N is equal to 80, unless otherwise noted. In 
our numerical analysis, we limit the particle number on 
each site and each species to two for filling v ^ 0.8 and 
four otherwise. Once we obtain the ground state, we cal- 
culate the energy, density distributions, correlation func- 
tions, and the structure factor to identify the quasi- long 
range order and other properties of the ground state. 

For example, to determine whether a SF, PSF, or 
CFSF is present, we study the decay behavior of the cor- 
relation functions, G{x), Ra{x), and Rs{x), defined in 
Eqs. El [HI and [H, respectively. If both Ra and Rs 
decay algebraically, the system is in a single-particle su- 
perfiuid (SF) state. If both are exponential, the system 
is in a Mott insulator(MI) state. If Rs or Ra decays 
algebraically, the system is in the PSF or CFSF state, 
respectively. These relationships are summarized in Ta- 
ble H 

In Fig. [5lja) and (b), we show the decay behavior of 
the correlation functions in the PSF and CFSF phase, 
respectively. As the Hamiltonian is discrete, the cor- 
relation functions are calculated as discrete functions: 

and 



(&I/a,,), Rs{l,j) = {h\A 



{b\,MM.3b\.i)- For the PSF phase, RA{i,j) 



decays exponentially, while Rs decays algebraically. It is 
also worthwhile to notice that the single-particle Green's 
function decays exponentially, implying the absence of 
single-particle superfiuidity. For the CFSF phase, Ra de- 
cays algebraically while Rs decays exponentially. Single- 
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Figure 5: The correlation functions Ra, Rs, and G on a log- 
arithmic scale as a function of distance \i — j\. The index i 
is 40, the center of the 80 lattice sites. The squares are the 
numerical data. The blue lines are exponential fits to the data 
and red dotted lines are algebraic fits. Note that the scale of 
the vertical axis of the graphs differs by orders of magnitude. 
In (a), we show an example for the paired superfiuid phase at 
1^=0.3, t = 0.02f/, and U12 — — 0.16f/. Ra decays exponen- 
tially and Rs decays algebraically. The single-particle corre- 
lation function decays exponentially, implying the absence of 
single-particle superfiuidity. In (b), we show an example for 
the counterfiow superfiuid phase at = 0.5, t = 0.02(7, and 
U12 — 0.2(7. The anti-pair correlation function Ra decays 
algebraically, while the pair correlation function decays expo- 
nentially. Single-particle superfiuidity is again absent. The 
algebraic fits deviate from the data around \i — j\ ~ 40, due 
to the boundary conditions of our numerical calculations. 



particle superfluidity is again absent. 

Behavior of Ks and Ka ■ We study the decay behav- 
ior of Rs and Ra in more detail. Using the flt function, 
c- \i — where c and a are the fitting parameters, 

we obtain the power-law exponent a and, hence, the Lut- 
tinger parameters Ks and Ka based on Eqs.fTTIandfTSl 
In Fig. [6l[a) , we show these Ks and Ka as a function of 
U12, for non-commensurate filling. A Luttinger parame- 
ter is formally set to zero when its correlation function 
decays exponentially. 

For U12 < — 0.06C/, Ra decays exponentially, while 
for U12 > — 0.06C/, Ra decays algebraically, and Ka- 
increases as U12 increases. The system undergoes a PSF 
to SF transition at U12 = — 0.06C/. On the other hand, 
Ks decreases monotonically for U12 > —0.6(7. For U12 < 
—0.6(7 the numerics failed to converge to a homogeneous 



state. This indicates that the system collapses, and we 
therefore cannot extract a Luttinger liquid parameter. 
We can observe charge density wave (CDW) order for a 
range of (7i2/(7 in Fig.[6l According to Eq. [161 this order 
exists when Ks + Ka < 2. In fact, it co-exists with the 
SF, PSF or CFSF order. At half-fiUing, Ks will go to 
zero at a critical, positive value of U12. This indicates 
the transition from the SF to CFSF phase. 

Finite-size effect: The behavior of Ka/s stated above 
is affected by the size of the system. Finite size effects 
can 'smooth out' a sudden change in Ka/s at the phase 
transition. This effect can be estimated from the RG 
flow calculation by integrating Eqs. [8] and [9] out to a fl- 
nite value I rather than to infinity. In Fig. (Hb) , we show 
an example of a finite-^ RG calculation in the vicinity of 
the PSF-to-SF transition. We see that as / increases, Ka 
dramatically changes for the attractive U12. In the limit 
of I 00, Ka becomes discontinuous and 'jumps' from 
to 1 at U12 » -0.01(7. This is where the PSF-to-SF tran- 
sition occurs. This transition is a Berezinskii-Kosterlitz- 
Thouless transition [l3,[2§|. In order to compare the RG 
result with our TEBD result, we associate the system 
size N with the flow parameter I, based on the relation 
in Eq. [TOl The cut-off tq is the lattice constant and 
Tq = Nql- For TV = 80 we have I = 4.4 and we find that 
the RG and TEBD are in good agreement. The regime 
between U12/U « —0.06 and —0.01 is a cross-over regime 
due to the finite size of the system. 

Collapse and phase separation: For large |(7i2|, the 
system approaches collapse or phase separation. Accord- 
ing to Tomonaga-Luttinger liquid theory, Ks ^ 00 as 
the system approaches collapse and Ka ^ 00 as the sys- 
tem approaches phase separation. As seen in Fig. El we 
indeed find such a tendency in our TEBD calculations. 
For U12 > 0.8(7 (not shown), Ka increases rapidly to 
values around 10, indicating a possible phase separation. 
For (7i2 < —0.6(7, due to the slow decay of the corre- 
lation function Rs and the finite-size of our system, we 
are unable to extract an accurate Ks from the numerical 
result. On the other hand, we observe a peaked density 
distribution for U12 < —0.6(7, indicating a collapse. In 
the phase separation regime, G(x) has algebraic decay 
except for ly — 0.5 or 1, where it has exponential de- 
cay. An algebraic decay implies two spatially-separated 
single-species superfiuids while the exponential decay im- 
plies two spatially-separated Mott insulators. (30|. 



A. Phase diagram 

We study the phase diagram as a function of filling j/ 
and parameters of the Hubbard Hamiltonian. Assuming 
a positive (7, the system can be fully characterized in 
terms of i^, t/(7, and U12/U. Our results are shown in 
Fig.|7]for a fixed hopping parameter and in Fig.|8]for half 
filling. 
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(a) 




Figure 6: (a) Ks and Ka as a function of U12 as extracted 
from the fit of the correlation functions, Rs and Ra- The 
filling V is 0.7 and t/U is 0.02. Around U12/U ^ -0.06, the 
anti-pair correlation function changes from algebraic to ex- 
ponential decay. This corresponds to the transition from the 
PSF to SF phase. When Ra decays exponentially, Ka is for- 
mally set to zero. For Ka + Ks < 2, the system has CDW 
order. Error bars are one standard deviation uncertainties 
obtained from the power-law fit to the numerical data, (b) 
A comparison of Ka obtained from our RG and TEBD cal- 
culations. The red square connected by lines are the TEBD 
results while all other lines are determined from the RG fiow 
with fiow parameter / = 3, 4, 7, and 10, where I is defined in 
Eq. [To] The error bars are as in panel (a). The PSF-to-SF 
transition obtained from TEBD is around U12/U ~ —0.06, 
while the RG calculation shows that for I = 10, the transition 
occurs near U12/U ~ —0.01. We interpret the regime between 
U12/U ~ —0.06 and —0.01 the cross-over region. 



1. Phase diagram at a fixed hopping parameter 

In Fig. [7] we show the phase diagram for filUng frac- 
tions between and 1 and the interaction U12/U between 
-1.1 and 1.1. The symbols correspond to numerical data 
points at which the phases have been characterized. Dif- 
ferent markers represent the different orders. The orders 
are determined from the decay behavior of the three cor- 
relation functions Ra, Rs, and G. 

For weak attractive inter-species interaction, —0.06 < 
U12/U < 0, the system is in a SF state. As U12 grows 
more attractive, paired superfluidity (PSF) occurs. The 
critical U12 is largest, — 0.08J7, at half-filling and grad- 



ually decreases away from half-filling. This phase bound- 
ary differs from that predicted by our RG calculation (Eq. 
[20l) . plotted as the dotted fine in Fig. [71 This discrepancy 
is the result of the finite-size effect discussed in Fig. [61(b). 
In the SF to PSF cross-over regime, charge density wave 
(CDW) order can coexist. According to the phase dia- 
gram Fig. [31 for attractive interaction, CDW order can 
co-exist only with PSF order. In our numerical work, we 
observed the CDW order slightly outside the numerical 
phase boundary of PSF but within the RG phase bound- 
ary of PSF. The sub-regime where CDW and PSF co- 
exist ends when U12/U < —0.4. When the inter-species 
attraction is comparable to the intra-species repulsion, 
U12 < —U, the system collapses (CL) and no long-range 
order is present. 

For repulsive inter-species interaction and U12 < U, 
the system is in a SF state for all non-commensurate 
fillings. Within the SF regime, there is a smaller param- 
eter region where CDW order coexist with the SF order. 
This subregime is a quasi-supersoHd regime. The bound- 
ary between a normal superfiuid and a quasi-supersolid 
is estimated by RG calculation in Ref. At half- 

fiUing, counterflow superfluidity (CFSF) occurs when 
0.08 < U12/U < 1. Within the CFSF regime, the 
CDW order can coexist, forming a quasi-supersolid of 
anti-pairs. It also worthwhile to point out that at half- 
flUing, CDW order only exists within the PSF and CFSF 
regimes. 

At unit fllling, our numerical results do not show ev- 
idence of PSF or CFSF for any U12. We flnd a Mott 
insulator (MI) state for \Ui2\ < U. 



2. Phase diagram at half-filling 

In Fig. [SI we show the phase diagram at half filling 
as a function of U12/U and t/U. From this diagram, 
we find that the border between PSF and SF and the 
border between PSF and CL approach each other as t 
increases. Similarly, the border between the CFSF and 
SF and the border between CFSF and PS approach each 
other. In fact, the PSF and CFSF phases end around 
t ~ 0.16C/. Within the PSF and CFSF regimes, CDW 
order can co-exist. In the phase separated regime, the 
separated single-species ensembles form two individual 
Mott insulating states for t < 0.1411 and two individual 
SF states for t > O.UU. 

We can compare this phase diagram with the half- 
filling phase diagram in Fig. [H obtained from Tomonaga- 
Luttinger liquid theory. Especially, we can compare 
the location of the phase boundary between SF and 
PSF(CFSF). To do so, we plot the RG phase bound- 
aries, described by Eqs. [201 and [211 onto our phase dia- 
gram. The area near the two boundaries is interpreted 
as the cross-over regime where finite-size effects modify 
the phase boundary. 
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Figure 7: Phase diagram for a homogeneous system with 80 sites and the hopping parameter t — 0.02U as a function of filling 
V and inter-species interaction U12/U . The horizontal axis shows three disconnected regions in U12/U . The solid lines are the 
estimated phase boundaries based on the TEBD results and the dotted line is the PSF-to-SF phase boundary predicted by our 
RG calculation (see Eg. I20p . For attractive interaction U12 ^ —0.0617, the system forms a paired-superfiuid (PSF). The state 
collapses(CL) for U12 < —0.7(7. For U12 > —0.06 and U12 < U the system shows single-particle superfluidity (SF). The system 
phase-separates (PS) for ?7i2 > 1 and forms two single-particle superfluids (SF). Open circles are the points where Ks + Ka < 2 
and charge density wave (CDW) order coexists with a superfluid phase (SF,PSF, or CFSF). At half and unit fllling there exist 
special phases. For repulsive interaction U12 > 0.08)7 and half-filling, the system forms a counterfiow superfiuid (CFSF). For 
unit filling, we find a Mott-Insulator (MI) phase for interactions \Ui2\ < U. Finally, in the PS region at half- and unit-filling, 
the system forms two individual MI states. 



B. Realization and detection 

Having established the phase diagram for the homo- 
geneous system, we now discuss how to reahze and de- 
tect the PSF and CFSF phases. First, we need to 
modify the Hubbard Hamiltonian in Eq. [1] because in 
any ultra-cold atom experiment an additional trapping 
potential is present. We add a harmonic potential, 
^{j ~ jcY(ni^j + "2,j), where j is the site index and 
jc is the index at the center of the system. The TEBD 
method is used to find the ground state. We consider a 
system of 80 lattice sites and adjust the total number of 
particles and the trap frequency so that the number of 
particles is negligible at the edge of the lattice. 

We again determine the orders of the system by study- 
ing the correlation functions in Table HI We find that, in 
spite of the presence of the trap, the correlation functions 
still show exponential or algebraic scaling away from the 
edge of the lattice. In fact, a correlation function can 
have different decay behavior in different parts of the 
trap. We also find that SF, PSF, and CFSF still exist. 
The remainder of this article focusses on experimental 
signatures that distinguish between these orders by cal- 



culating the density distibution, the time-of-fiight image 
after an expansion, or the structure factor for Bragg spec- 
troscopy. 

Density distribution: We find that in a trapped system 
PSF and CFSF can only exist when the density distri- 
bution satisfies certain conditions. For PSF, the density 
of each species at the center of the trap, ^center, must 
be less than one atom per site or equivalently per lattice 
constant a^. (The density is largest at the center.) For 
CFSF, ncontcr must satisfy nccntcrCiL = 1/2. Once such 
conditions are satisfied, the critical value of U12 for PSF 
and CFSF is close to the one for a homogeneous system 
(See Figs. [7] and [HI . 

In Fig.[9lja) we show density distributions for three at- 
tractive interactions U12 and a hopping parameter equal 
to the one used for Fig. [71 For all attractive interac- 
tions, the density distributions of each species are the 
same. For more attractive inter-species interaction, the 
density distribution concentrates near the center of the 
trap. There is no discontinuous change in the density 
distribution when the system goes from SF to PSF. 

In Fig.[9Ub) we show the density distribution for f/i2 = 
Q.2U . In this case in the center of the trap, where the 
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Figure 8: Phase diagram at half-filling as a function of Uu/U and t/U . The solid lines are estimated phase boundaries from 
the TEBD calculation and the dotted lines are the phase boundaries predicted by the RG calculation (see Egs. I20l and l24p . For 
large repulsive interaction, the system phase separates (PS) and for large attractive interaction, the system collapses (CL). For 
moderate interactions and for t/U < 0.2, the system shows paired superfluidity (PSF) on the attractive side and counterflow 
superfluidity on the repulsive side. Both PSF and CFSF can coexist with charge density wave (CDW) order when t < Q.IU . 



density distribution is constant or has a "plateau", the 
system is in a CFSF state. The "plateau" is at half-filling 
consistent with predictions from a local density approx- 
imation and noting that in Fig. [7] CFSF only occurs at 
V = 1/2. Towards the edge, where the density is decreas- 
ing sharply, it is in a SF state. The plateau impHes that 
the system is incompressible in the center. 

Time of flight measurement: A widely used measure- 
ment technique in the field of ultra-cold atoms is mea- 
suring the density of atoms after a time-of-flight (TOF) 
expansion. The ID optical lattice potential and the har- 
monic trap are abruptly turned off at time T — Q and 
the atoms expand freely afterwards. We calculate the 
density at time T, according to 

na{x,T) = {ci{x,T)ca{x,T)) (31) 

with a = 1,2. The operators Ca{x,T) are related to the 
lattice operator ba,j according to 

N 

ha{x,T) = Y,w{x-rj,T)ba.j, (32) 
i=i 

where w(x,T) = ^Jd/V2^A{T)^ exp(-a;V(4A(T)2)) de- 
scribes the free expansion from the initial Gaussian wave- 
function of an atom in a lattice site and A(T)^ = 
(P + iTh/{2m). The parameter d is the width of the 
initial Gaussian state and m is the atomic mass. The 



density distribution na{x,T) is then given by 

N 

na{x,T)^ ^ w*{x-rj^,T)w{x-rj^,T)G{ji,j2), 

j 1,32 = 1 

where G(ji, j2) is the single-particle Green's function. In 
Fig. [10] we show examples of TOF expansions of PSF, 
CFSF, and SF order. For the SF phase, we find a strongly 
peaked interference pattern, refiecting the single-particle 
quasi-long range order. For both PSF and CFSF phases, 
the TOF density shows a broad Lorentzian distribution, 
which is due to the exponential decay of the single- 
particle Green's function. 

Feshbach ramp: In order to detect the superfluidity 
of pairs, we consider applying a Feshbach ramp to pair- 
wise project the atoms onto molecules formed by one 
atom from each species, which is similar to detection 
of fermionic pairs in the BCS regime 01 • In those ex- 
periments, a fast ramp across a Feshbach resonance was 
used, followed by a time-of-flight expansion. The den- 
sity distribution of the molecules showed the superfluid- 
ity of fermionic pairs. We propose a similar detection for 
bosonic pairs in PSF. 

To give a simple estimate of a TOF image after a Fesh- 
bach ramp, we imagine that bosons of different species on 
the same lattice site are converted into molecules. This 
leads to the replacement ^i,j&2,j Mj, where Mj is the 
molecule annihilation operator. A TOF density of the 
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Figure 9: Density distribution of a trapped system for t = 
0.02(7. (a) Attractive interaction C/12. The trap frequency is 
Q = Ix 10~^ U and the number of atoms is 20 for each species. 
For attractive interactions, the density distributions of the 
two species are identical. For U12 — —0.01(7 (curve I) the 
system is superfluid. For U12 = —0.11(7 (curve II) and U12 = 
— 0.21U (curve III), the system is in the paired superfluid 
(PSF) state. As U12 becomes more negative the distribution 
gradually shrinks in size, (b) Repulsive interaction U12 = 
0.2(7 with Q = 8 X lO'^U and 30 atoms of each species. The 
red and green curves correspond to the species, respectively. 
The density distribution has a 'plateau' with half-flUing in the 
center of the trap. The system is in a counter-flow superfluid 
(CFSF) state. The two species have weak interlocked density 
modulations around half filling. 




Figure 10: Density distribution after a time-of-flight expan- 
sion. We assume ^'^Rb atoms and use an expansion time of 
0.03s. The hopping energy is t = 0.02(7. Panel (a): For 
attractive interaction U12, we show the TOF expansion of a 
SF state at U12 = -0.01(7 (red line) and of a PSF state at 
U12 = —0.21(7 (green line). The two curves correspond to 
the expansion of the densities shown as curve I and III in 
Fig. M,&) The trap frequency is Q = 1 x 10~^U. Panel (b): 
For repulsive interaction, we show a TOF expansion of a SF 
state at U12 = 0.01(7 and of a CFSF state at U12 = 0.21(7. 
The trap frequency is f2 = 8 x 10~^I7. 



molecules at position x and time T is given by 

JV 

nM{x,T)= ^ w*{x - rj^,T)w{x - rj.,_,T)Rs{ji,i2)- 

31,32 = 1 

(33) 

In the expanding wave function ■w{x,T), the mass m is 
replaced by the mass of the molecule. We assume the 



same initial vifidth d. In a more realistic estimate, the 
conversion efRciency to molecules would not be 100%, 
but approximately given by the square of the overlap of 
the molecular wave function and the single-atom wave 
functions. This leads to a reduced signal. The spatial 
dependence, however, remains the same. In Fig. [iT], we 
see an example of the density of molecules after TOF 
and, for comparison, the atomic density after TOF for the 
PSF state. The strongly peaked molecular distribution 
indicates the quasi-condensate of the bosonic pairs. The 
single-atom density is a broad Lorentzian distribution, 
indicating the absence of single-particle SF. 

Bragg spectroscopy: To detect the presence of CDW 
order, one can use Brag spectroscopy jsil,!!^!- The quan- 
tity that is measured in those experiments is either the 
dynamic or static structure factor. Here we calculate the 
static structure factor Sa{k) for species a = 1,2. It is 
defined as 



Saik) 



N ^ 

31 ,32 



-ikaL{3l-32) 



{{na{jl)na{32)) 



{na{jl)){na{j2)))- 



(34) 



For wavevectors k near twice the "Fermi wavevector" 
kp, the structure factor S'(fc) ~ ||fc| - 2fcFr""'='""'' with 
acDW = 2 — Ks — Ka In our system, Ks + Ka 

is always larger than 1 and, thus, 1 — acDW is positive. 
Consequently, the structure factor does not diverge. In 
the CDW regime with Ks + Ka < 2 the power l — acDW, 
however, is less than one. This gives S{k) cusps at ±2fcF 
when CDW quasi-long range order is present. In Fig. [12] 
we show examples of S{k) for a case with and without 
CDW. 

Bragg Spectroscopy preceded by a tt/2 pulse: To de- 
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Figure 11: Density distribution of molecules after time-of- 
flight expansion of state III in Fig. [9Ua). The expansion time 
is 0.03s. We assume two hyperfine states of Rb. These are 
converted into Feshbach molecules at T = via a fast ramp 
across a resonance. We assume a complete conversion. The 
strongly peaked interference pattern of molecules indicates 
the presence of a quasi-condensate of pairs. For comparison, 
we also show the TOF expansion of atoms in the PSF phase 
for the same parameters. The broad Lorentzian distribution 
demonstrates the absence of single-particle SF. 
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Figure 12: Structure factor at filling v — 0.3. For U\2 — 
— O.OIC/ the system is in the SF regime (dashed line) and for 
Vvi = —0.07(7 the system is in the PSF regime (continuous 
line). Cusps at |fc| = 2-kv only occur for (7i2 = —0.0717 indi- 
cating the coexistence of CDW with PSF order. 



tect CFSF order, we propose the following detection 
method. It applies to the case that the mixture is 
composed of atoms in different internal states rather 
than different atomic species. First, we apply a 7r/2 
pulse, which transfers the atoms into the superposi- 
tions ^1/2,1 b±^i — ± hi.i)l\f2. We then mea- 
sure the structure factor, which now corresponds to the 
Fourier transform of the density correlations i?n±(j, j) — 
{n±^in±j) — {n±^i){n±j). In terms of the original foi/2.i 
operators these density correlations are given by 



both Tomonaga-Luttinger liquid theory and the time- 
evolving block decimation method. We first discussed 
the zero-temperature phase diagram in a homogeneous 
system at different filling fractions and different param- 
eter regimes. We have shown that ID Bose mixtures in 
an optical lattice can have quasi-long range orders that 
include superfluid, paired superfluid (PSF), counterflow 
superfiuid (CFSF), and Mott insulator. We also found 
that each type of superfiuid order can coexist with charge 
density wave (CDW) order and that in both PSF and 
CFSF phases single particle superfiuidity (SF) is absent. 



In addition, we discussed ways of realizing and de- 
tecting these phases experimentally. We propose using 
a Feshbach ramp to probe the momentum distribution of 
pairs in the PSF, which shows signatures of the quasi- 
condensate of pairs. To detect the CFSF for a mixture 
composed of two atomic hyperfine states, we propose to 
measure the static structure factor by using Bragg spec- 
troscopy preceded by a 7r/2 pulse. A sharp peak in the 
structure factor was shown to be dominated by the con- 
tribution from the momentum distribution of anti-pairs 
in the CFSF phase. Finally, we suggest to detect CDW 
order with Bragg spectroscopy. 
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({fii^i + n2,i){ni,j + n2.j)) 



(35) 



The last term in the above equation is the correlation 
function i?a(«jj) of the order parameter of CFSF, bijblj. 
In Fig. [I3l we show the structure factor S+{k), the 
Fourier transform of Eq. [35l as well as the Fourier trans- 
form of Ra{i,j). Both S+{k) and the Fourier transform 
of Ra{i,j) have a cusp around fc = 0. The cusp is due to 
the long-range correlations of the anti-pairs in the CFSF. 
The two functions are nearly identical near k = 0, indi- 
cating that the momentum distibution of anti-pairs can 
be measured by determining the structure factor S+{k) . 
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V. SUMMARY 

We have studied ground state properties of one- 
dimensional Bose mixtures in an optical lattice using 



Figure 13: Structure factor S+{k) (blue line) after applying a 
7r/2 pulse in the CFSF phase. The quasi-condensate of anti- 
pairs generates an algebraic peak at fc = 0. The cusp also 
appear in the Fourier transform of the anti-pair correlation 
function Ra{i,j) = {b\ ib2,ib2 jbi,j){ied dashed line). 
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Appendix A: TEBD METHOD FOR 
TWO-SPECIES MANY-BODY SYSTEMS 

In this appendix, we briefly review the time-evolving 
block decimation (TEBD) method used in Sec. |EI 
and explain an efficient way to apply the TEBD to a 
two-species Bose-Hubbard model. We use the number- 
conserving version of the TEBD method [lB |. 

The TEBD determines the ground state via an imag- 
inary time evolution for one-dimensional (ID) quantum 
lattice systems. In this method the Hilbert space H is 
decomposed as 

n = ®fijii. (Ai) 

Here, I refers to the /th lattice site, M is the number of 
sites, and H/ is the local Hilbert space at site I with local 
dimension d, independent of /. Any state \^) in H is 
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represented as 

d 

l*>- c,i,,„...,,„|ji)|j2>---|jAf>. (A2) 

il J2,--- jAf = l 

In the TEBD algorithm, coefficients Cj-^j^^,,,^^^ are de- 
composed as 

XI X2 Xm- 1 

01 — 102 — 1 Oj\-/ — 1 — 1 

■V [M-2] p [M- 1] jAf - 1 \ [M- 1] p lM]jM (Ao\ 

'^'*aAf-2 0Ar_20Af_l 'OAJ-I OAf-l ' V^^V 

The variables Aq| and x/ are the Schmidt coefficients and 
rank of the Schmidt decomposition of l^*) with respect to 
the bipartite splitting of the system into [1, . . . , ^ — 1, ^] : 

l*)-EASl$LV-''-'''')l*lir''+''-'''')- (A4) 

oi=l 

We take X^^ > a)^' for all a < (3. In one dimension, 
the rank xi at the center of the system must be of the 
order d*^/^ in order to express arbitrary states. However, 
since it is empirically known that the Schmidt coefficients 
Ao decrease rapidly with index a for the ground and low- 
lying excited states, we set xi to a relatively small number 
X for all I. 

To efficiently simulate the two-species Bose-Hubbard 
model (Eq. [T]in the main text), we map it onto the one- 
species Hamiltonian 

2N-2 

H = -t ibjbi+2 + h.c.) + Ui2 Y '^I'^i+i 

1 = 1 odd J 

u 

1=1 

where N is the number of sites in the original two-species 
Hamiltonian. In this one-species Hamiltonian, there are 
2N sites, each of which is indexed by I. The odd sites 
I correspond to species 1 and the even sites to species 
2. Hopping between neighboring sites —tb'l fia,i+i in 
Eq. [T] is mapped onto a next-nearest-neighbor hop- 
ping —tblbi+2 in Eq. IA5I Similarly, the inter-species 



onsite-interaction Ui2ni^in2,i is mapped onto the nearest- 
neighbor interaction Uunini+i. This type of mapping 
has been successfully applied to treat the two-legged 
Bose-Hubbard model [l7| . 

We map the two-species Bose-Hubbard Hamiltonian 
Eq.[T]onto the one-species Hamiltonian because it reduces 
computational cost dramatically. This cost in TEBD jTl| 
scales as Md^x'^- For the two-species system with N sites 
we must define a dimension of the local Hilbert space for 
each species, say D. Hence, at each site there are 
basis functions and the cost scales as ND^ . On the other 
hand, for the mapped Hamiltonian with 2N sites and a 
local dimension D the cost only scales as 2ND^. In our 
calculation, we set d = 3 for the ffiling factor ly < 0.8 and 
d = 5 for ly = 0.9, 1. In this case, the mapping makes the 
computation five to ten times faster. 

Imaginary time evolution of any state to the ground 
state is given by repeated application of e~^^^ on \^), 
where i5 is a small imaginary time step. To apply this 
operator we first split the Hamiltonian into three parts 
as = i/int + H°^^ + HZ^, where 



N 

Hint = Y [^12 
m— 1 

+ Un2min2m - 1)] , (A6) 
^hop = II (6L-ife2,„+i + &LW2+h.c.), 

odd rn 

H^lp = E (4n-lfe2™+l+6LW2+h.C.). 

even m 

Subsequently, we use the second-order Suzuki- Trotter ex- 
pansion to decompose e~'^^^ as 

e-'"^ = g-«-ff.nt5/2g-»//Sop'5/2g-»/fwop"«e-''-f^hop'5/2 

xe-''^'»'*/2 + 0(j3), (A7) 

Each of the operators e"'^'-"^/^^ ^--^H^ap^/^^ and 
g-'-H^hop "5 (-an be decomposed into a product of two-site 
operators, which can be efficiently applied to the ma- 
trix product state |^) (l3 . ITsl . [la |. We use swapping 
techniques to apply the next-nearest-neighbor operators 
g~^^^.tp^5/2 ande-^^-p"* [11,0. 



